!############################# Change Log ##################################
! 5.0.0
!
!###########################################################################
!  Copyright (C)  1990, 1995, 1999, 2000, 2003 - All Rights Reserved
!  Regional Atmospheric Modeling System - RAMS
!###########################################################################

subroutine diffuse()

  ! +-----------------------------------------------------------------+
  ! \	this routine is the subdriver to compute tendencies due to    \
  ! \	  subgrid-scale turbulence.				      \
  ! +-----------------------------------------------------------------+

  use mem_tend, only:    &
       tend                   ! %tket, %epst, %ut, %vt, %wt

  use mem_basic, only:   &
       basic_g                !  %up(IN), %vp(IN)

  use var_tables, only:  &
       num_scalar,       &    ! INTENT(IN)
       scalar_tab             ! %var_p, %var_t

  use mem_turb, only:    &
       idiffk,           &    !INTENT(IN)
       turb_g,           &    ! %tkep, %hkm, %vkh
       xkhkm

  use mem_grid, only:     &
       if_adap,           &  !           INTENT(IN)
       jdim,              &  !           INTENT(IN)
       ngrid,             &  !           INTENT(IN)
       grid_g,            &  ! %rtgt     INTENT(IN)
       dzm,               &  !           INTENT(IN)
       dzt,               &  !           INTENT(IN)
       npatch,            &  !           INTENT(IN)
       nstbot,            &  !           INTENT(IN)
       nscl,              &  !           INTENT(IN)
       naddsc                !           INTENT(IN)

  use mem_leaf, only:     &
       leaf_g                  !INTENT(IN)

  use mem_micro, only:    &
       micro_g                 !%rcp

  use mem_scratch, only:  &
       scratch,           &    ! %vt3da, %vt3db, %vt3dc, %vt3dd, %vt3de,
                               ! %vt3df, %vt3dg, %vt3dh, %vt3di, %vt3dj,
                               ! %vt3dn, %scr2
       vctr34                  !

  use node_mod, only:     &
       mxp,               &  !INTENT(IN)
       myp,               &  !INTENT(IN)
       mzp,               &  !INTENT(IN)
       ia,                &  !INTENT(IN)
       iz,                &  !INTENT(IN)
       ja,                &  !INTENT(IN)
       jz,                &  !INTENT(IN)
       ia_1,              &  !INTENT(IN)
       ja_1,              &  !INTENT(IN)
       iz1,               &  !INTENT(IN)
       jz1,               &  !INTENT(IN)
       ibcon,             &  !INTENT(IN)
       mynum,             &  !INTENT(IN)
       mi0,               &  !INTENT(IN)
       mj0,               &  !INTENT(IN)
       mmyp,              &  !INTENT(IN)
       mmxp,              &  !INTENT(IN)
       mmzp,              &  !INTENT(IN)
       izu,               &  !INTENT(IN)
       jzv,               &  !INTENT(IN)
       ia1,               &  !INTENT(IN)
       ja1,               &  !INTENT(IN)
       iz_1,              &  !INTENT(IN)
       jz_1                  !INTENT(IN)

  use ke_coms, only:      &
       alf_eps,           &
       alf_tke

  use micphys, only:      &
       level

  use mem_turb_scalar, only:   &
       turb_s

  ! CATT
  use mem_grell,  only:        &
       grell_g_sh

  use catt_start, only:        &
       CATT           ! intent(in)
  
  use mem_cuparm, only : &
       nnqparm

!--(DMK-CCATT)---------------------------------------------------------
  !srf-chem
  use chem1_list
  use mem_chem1, only: &
       chem1_g, &
       CHEMISTRY, &
       NSPECIES_TRANSPORTED
  !srf-chem-end
  !ml/srf- for stilt/new turb scheme
  use mem_stilt
!--(DMK-CCATT-END)-----------------------------------------------------

  implicit none

  !local variables:
  integer :: mxyzp,ind,n
  real :: s1,s2,s3
  real, pointer :: scalarp, scalart,vkh_p,hkh_p
  integer :: i,j,k,ksf

  !###########################################################################
  ! CATT
  ! srf - fev-2003 - Large Scale Forcing for shallow and deep cumulus

!--(DMK-CCATT)---------------------------------------------------------
  real, pointer :: lsfcupar_p
!--(DMK-original)------------------------------------------------------
!  real, pointer :: lsfcupar_p(:,:,:)
!--(DMK-CCATT-END)-----------------------------------------------------

  !###########################################################################

  ! Interface to diffsclr needed because of a pointer argumenter
!!$  interface
!!$     subroutine diffsclr(m1,m2,m3,ia,iz,ja,jz,jd  &
!!$          ,ia_1,ja_1,ia1,ja1,iz_1,jz_1,iz1,jz1,n,ksf  &
!!$          ,scp,sct,vt3da,vt3db,vt3df,vt3dg  &
!!$          ,vt3dj,vt3dk,vt3do,vt3dc,dn03i,vt3dl,vt3dm,vt2db,rtgt,sfcflx  &
!!$          ,dn0,vkkh,hkkh, &
!!$          ! CATT - srf - Large Scale Forcing for GRELL CUPAR
!!$          ! LFR      iscalar,lsfcupar)
!!$          lsfcupar,nsc)
!!$       integer, intent(IN) :: n,ksf,nsc
!!$       integer, INTENT(IN) :: m1    &
!!$            , m2    &
!!$            , m3    &
!!$            , ia    &
!!$            , iz    &
!!$            , ja    &
!!$            , jz    &
!!$            , jd    &
!!$            , ia_1  &
!!$            , ja_1  &
!!$            , ia1   &
!!$            , ja1   &
!!$            , iz_1  &
!!$            , jz_1  &
!!$            , iz1   &
!!$            , jz1
!!$       real, dimension(m1,m2,m3), INTENT(IN) :: scp     &
!!$            , vkkh    &
!!$            , hkkh    &
!!$            , dn0
!!$       real, dimension(m1,m2,m3), INTENT(INOUT) :: sct    &
!!$            , vt3da  &
!!$            , vt3db  &
!!$            , vt3dc  &
!!$            , vt3df  &
!!$            , vt3dg  &
!!$            , vt3dj  &
!!$            , vt3dk  &
!!$            , vt3do  &
!!$            , dn03i  &
!!$            , vt3dl  &
!!$            , vt3dm
!!$       real, dimension(m2,m3), INTENT(IN)    :: sfcflx,rtgt
!!$       real, dimension(m2,m3), INTENT(INOUT) :: vt2db
!!$       ! CATT
!!$       !srf- Large Scale Forcing for GRELL CUPAR
!!$       !integer ::  iscalar
!!$       real, pointer :: lsfcupar(:,:,:) ! dimension(m1,m2,m3)
!!$       !srf- ------------------------------------
!!$     end subroutine diffsclr
!!$  end interface

  ! CATT
  ! Nullifing pointer to Large Scale Forcing for GRELL CUPAR - Not used
  nullify(lsfcupar_p)
  !

  mxyzp = mxp * myp * mzp

  scratch%vt3dg = 0. !LFR - Avoiding overflow in strain

  if (if_adap == 0) then

     call strain(mzp,mxp,myp,ia,iz,ja,jz                       &
          ,ia_1,ja_1,iz1,jz1,jdim                                &
          ,basic_g(ngrid)%up (:,:,:) ,basic_g(ngrid)%vp (:,:,:)  &
          ,basic_g(ngrid)%wp (:,:,:) ,scratch%vt3da     (:)      &
          ,scratch%vt3db     (:)     ,scratch%vt3dc     (:)      &
          ,scratch%vt3dd     (:)     ,scratch%vt3de     (:)      &
          ,scratch%vt3df     (:)     ,scratch%vt3dg     (:)      &
          ,scratch%vt3dh     (:)     ,scratch%vt3di     (:)      &
          ,scratch%vt3dn     (:)     ,scratch%scr2      (:)      &
          ,idiffk(ngrid))

  else

     call strain_adap(mzp,mxp,myp,ia,iz,ja,jz                  &
          ,ia_1,ja_1,iz1,jz1,jdim                                &
          ,grid_g(ngrid)%lpu (:,:)   ,grid_g(ngrid)%lpv (:,:)    &
          ,grid_g(ngrid)%lpw (:,:)   ,basic_g(ngrid)%up (:,:,:)  &
          ,basic_g(ngrid)%vp (:,:,:) ,basic_g(ngrid)%wp (:,:,:)  &
          ,scratch%vt3da     (:)     ,scratch%vt3db     (:)      &
          ,scratch%vt3dc     (:)     ,scratch%vt3dd     (:)      &
          ,scratch%vt3de     (:)     ,scratch%vt3df     (:)      &
          ,scratch%vt3dg     (:)     ,scratch%vt3dh     (:)      &
          ,scratch%vt3di     (:)     ,scratch%vt3dn     (:)      &
          ,scratch%scr2      (:)     ,idiffk(ngrid)              &
          ,grid_g(ngrid)%dxm (:,:)   ,grid_g(ngrid)%dxt (:,:)    &
          ,grid_g(ngrid)%dxu (:,:)   ,grid_g(ngrid)%dxv (:,:)    &
          ,grid_g(ngrid)%dym (:,:)   ,grid_g(ngrid)%dyt (:,:)    &
          ,grid_g(ngrid)%dyu (:,:)   ,grid_g(ngrid)%dyv (:,:)    &
          ,dzm,dzt)

  endif

  if (level <= 1) call azero(mxyzp,scratch%vt3dp(:))
  if (level >= 2) call ae1  (mxyzp,scratch%vt3dp(:),micro_g(ngrid)%rcp(:,:,:))

!  call bruvais(mzp,mxp,myp,ia,iz,ja,jz 			 &
!	,basic_g(ngrid)%theta (:,:,:) ,basic_g(ngrid)%rtp (:,:,:)  &
!	,basic_g(ngrid)%rv    (:,:,:) ,scratch%vt3dp(:) 	   &
!	,basic_g(ngrid)%pp    (:,:,:) ,basic_g(ngrid)%pi0 (:,:,:)  &
!	,scratch%vt3dj        (:)     ,grid_g(ngrid)%rtgt (:,:)    &
!	,grid_g(ngrid)%lpw    (:,:)   )


!-srf 29/12/2008 adapted from OLAM 
  call bruvais_OLAM(mzp,mxp,myp,ia,iz,ja,jz			     &
       ,basic_g(ngrid)%theta (1,1,1) ,basic_g(ngrid)%rtp (1,1,1)  &
       ,basic_g(ngrid)%rv    (1,1,1) ,scratch%vt3dp(1)  	  &
       ,basic_g(ngrid)%pp    (1,1,1) ,basic_g(ngrid)%pi0 (1,1,1)  &
       ,scratch%vt3dj	     (1)     ,grid_g(ngrid)%rtgt (1,1)    &
       ,grid_g(ngrid)%lpw    (1,1)   )


!--(DMK-CCATT)---------------------------------------------------------
  !ml/srf- for new turn scheme
  if (idiffk(ngrid) <= 3 .or. idiffk(ngrid) == 7) then
!--(DMK-original)------------------------------------------------------
!  if (idiffk(ngrid) <= 3) then
!--(DMK-CCATT-END)-----------------------------------------------------

     call mxdefm(mzp,mxp,myp,ia,iz,ja,jz,ibcon,jdim            &
          ,scratch%vt3dh      (:)     ,scratch%vt3di      (:)    &
          ,scratch%vt3dj      (:)     ,scratch%vt3dk      (:)    &
          ,scratch%scr1       (:)     ,scratch%scr2       (:)    &
          ,basic_g(ngrid)%dn0 (:,:,:) ,grid_g(ngrid)%rtgt (:,:)  &
          ,grid_g(ngrid)%dxt  (:,:)   ,grid_g(ngrid)%dyt  (:,:)  &
          ,grid_g(ngrid)%lpw  (:,:)   ,mynum  )

     ! CATT
     if (CATT == 1) then
        !srf------
        !coef de difusao horizontal diferente (dum4 e ivt3dp) para tracers
        call mxdefm_tracer(mzp,mxp,myp,ia,iz,ja,jz  &
             ,ibcon,jdim,scratch%vt3dh(:),scratch%scr3(:) &
             ,basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%dxt(:,:),&
             grid_g(ngrid)%dyt(:,:),grid_g(ngrid)%lpw(:,:),mynum)

        !srf------
     endif

  endif

!--(DMK-CCATT)---------------------------------------------------------
  !ML -> Nananishi and Niino (2004) scheme based on Mellor-Yamada Level 2.5
  if (idiffk(ngrid) == 7) then
     call nakanishi(mzp, mxp, myp, ia, iz, ja, jz, jdim                       &
        ,turb_g(ngrid)%tkep       (1,1,1)   ,tend%tket                (1)     &
        ,scratch%vt3dd            (1)       ,scratch%vt3de            (1)     &
        ,scratch%vt3dh            (1)       ,scratch%vt3di            (1)     &
        ,scratch%vt3dj            (1)       ,scratch%scr1             (1)     &
        ,grid_g(ngrid)%rtgt       (1,1)     ,basic_g(ngrid)%theta     (1,1,1) &
        ,basic_g(ngrid)%rv        (1,1,1)   ,basic_g(ngrid)%dn0       (1,1,1) &
        ,basic_g(ngrid)%up        (1,1,1)   ,basic_g(ngrid)%vp        (1,1,1) &
        ,leaf_g(ngrid)%patch_rough(1,1,1)   ,leaf_g(ngrid)%tstar      (1,1,1) &
        ,leaf_g(ngrid)%ustar      (1,1,1)   ,leaf_g(ngrid)%patch_area (1,1,1) &
        ,turb_g(ngrid)%sflux_u    (1,1)     ,turb_g(ngrid)%sflux_v    (1,1)   &
        ,turb_g(ngrid)%sflux_t    (1,1)     ,grid_g(ngrid)%lpw        (1,1)   &
        ,stilt_g(ngrid)%pblhgt    (1,1)     ,stilt_g(ngrid)%lmo       (1,1)   &
        ,stilt_g(ngrid)%ltscale   (1,1,1)   ,stilt_g(ngrid)%sigw      (1,1,1))
  endif
  !ML
!--(DMK-CCATT-END)-----------------------------------------------------

  if (idiffk(ngrid) == 1) then
     call tkemy(mzp,mxp,myp,ia,iz,ja,jz,ibcon,jdim,mi0(ngrid),mj0(ngrid)  &
          ,turb_g(ngrid)%tkep   (:,:,:) ,tend%tket            (:)      &
          ,scratch%vt3dh        (:)     ,scratch%vt3di        (:)      &
          ,scratch%vt3dj        (:)     ,scratch%scr1         (:)      &
          ,grid_g(ngrid)%rtgt   (:,:)   ,basic_g(ngrid)%theta (:,:,:)  &
          ,basic_g(ngrid)%dn0   (:,:,:) ,basic_g(ngrid)%up    (:,:,:)  &
          ,basic_g(ngrid)%vp    (:,:,:) ,basic_g(ngrid)%wp    (:,:,:)  &
          ,turb_g(ngrid)%sflux_u(:,:)   ,turb_g(ngrid)%sflux_v(:,:)    &
          ,turb_g(ngrid)%sflux_w(:,:)   ,turb_g(ngrid)%sflux_t(:,:),vctr34 &
          ,grid_g(ngrid)%lpw    (:,:)   ,grid_g(ngrid)%lpu    (:,:)   &
          ,grid_g(ngrid)%lpv    (:,:))
  endif

  if (idiffk(ngrid) == 4) then
     call mxtked(mzp,mxp,myp,ia,iz,ja,jz  &
          ,ibcon,jdim  &
          ,turb_g(ngrid)%tkep   (:,:,:) ,tend%tket            (:)      &
          ,basic_g(ngrid)%up    (:,:,:) ,basic_g(ngrid)%vp    (:,:,:)  &
          ,basic_g(ngrid)%wp    (:,:,:) ,basic_g(ngrid)%rtp   (:,:,:)  &
          ,basic_g(ngrid)%rv    (:,:,:) ,basic_g(ngrid)%theta (:,:,:)  &
          ,scratch%vt3da        (:)     ,scratch%vt3dc        (:)      &
          ,scratch%vt3dh        (:)     ,scratch%vt3dj        (:)      &
          ,scratch%scr1         (:)     ,scratch%scr2         (:)      &
          ,turb_g(ngrid)%sflux_u(:,:)   ,turb_g(ngrid)%sflux_v(:,:)    &
          ,turb_g(ngrid)%sflux_w(:,:)   ,turb_g(ngrid)%sflux_t(:,:)    &
          ,grid_g(ngrid)%dxt    (:,:)   ,grid_g(ngrid)%rtgt   (:,:)    &
          ,grid_g(ngrid)%lpw    (:,:)   )
  endif

  !_STC............................................................
  !_STC Call to subroutine tkescl for E-l closure
  !_STC (S. Trini Castelli)
  !_STC............................................................
  if (idiffk(ngrid) == 5) then
     call tkescl(mzp,mxp,myp,npatch,ia,iz,ja,jz  &
          ,turb_g(ngrid)%tkep(:,:,:),tend%tket(:)  &
          ,turb_g(ngrid)%epsp(:,:,:),tend%epst(:)  &
          ,scratch%vt3da(:),scratch%vt3dc(:)  &
          ,scratch%vt3dh(:),scratch%vt3di(:)  &
          ,scratch%vt3dj(:),scratch%scr1(:)  &
          ,scratch%scr2(:) ,grid_g(ngrid)%rtgt(:,:)  &
          ,scratch%vt3dd(:),scratch%vt3de(:),grid_g(ngrid)%dxt(:,:)  &
          ,leaf_g(ngrid)%ustar(:,:,:),leaf_g(ngrid)%patch_area(:,:,:) &
          ,grid_g(ngrid)%lpw(:,:),basic_g(ngrid)%dn0(:,:,:)  )
  endif
  !_STC............................................................
  !_STC Call to subroutine tkeeps for E-eps closure
  !_STC (S. Trini Castelli)
  !_STC............................................................
  if (idiffk(ngrid) == 6) then
     call tkeeps(mzp,mxp,myp,npatch,ia,iz,ja,jz  &
          ,turb_g(ngrid)%tkep(:,:,:),tend%tket(:)  &
          ,turb_g(ngrid)%epsp(:,:,:),tend%epst(:)  &
          ,scratch%vt3da(:),scratch%vt3dc(:)  &
          ,scratch%vt3dh(:),scratch%vt3di(:)  &
          ,scratch%vt3dj(:),scratch%scr1(:)  &
          ,scratch%scr2(:) ,grid_g(ngrid)%rtgt(:,:)  &
          ,leaf_g(ngrid)%ustar(:,:,:),leaf_g(ngrid)%patch_area(:,:,:) &
          ,grid_g(ngrid)%lpw(:,:),basic_g(ngrid)%dn0(:,:,:)  )
  endif
  !_STC..................................................
  !_STC    Note: from subroutines TKESCL, TKEEPS :
  !_STC           VT3DI=Ke
  !_STC           SCR1=Km
  !_STC           VT3DH = Kh
  !_STC           SCR2 = SCR1 = Km
  !_STC..................................................
  !_STC............................................................

  call klbnd(mzp,mxp,myp,ibcon,jdim  &
       ,scratch%scr1 (:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:))
  call klbnd(mzp,mxp,myp,ibcon,jdim  &
       ,scratch%scr2 (:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:))
  call klbnd(mzp,mxp,myp,ibcon,jdim  &
       ,scratch%vt3dh(:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:))


  ! CATT
  if (CATT == 1) then
     !srf----
     !LFR call klbnd(mzp,mxp,myp,ibcon,jdim,a(ivt3dp),a(idn0))
     !call klbnd(mzp,mxp,myp,ibcon,jdim,scratch%vt3dp(:) &

!!$     call klbnd(mzp,mxp,myp,ibcon,jdim,scratch%vt3dp(:) &
!!$          ,basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:))

     call klbnd(mzp,mxp,myp,ibcon,jdim,scratch%scr3(:) &
          ,basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:))

     !srf----

  endif

  !_STC ....... boundary conditions even on Ke diffusion coefficient
  if(idiffk(ngrid) ==  5 .or. idiffk(ngrid) == 6) &
       call klbnd(mzp,mxp,myp,ibcon,jdim  &
       ,scratch%vt3di(:),basic_g(ngrid)%dn0(:,:,:),grid_g(ngrid)%lpw(:,:))

  !bob  swap new hkm, vkm, and vkh with past time level:  lagged K's have
  !bob  internal lateral boundary values from neighboring nodes

  ind = 0
  do j = 1,mmyp(ngrid)
     do i = 1,mmxp(ngrid)
        do k = 1,mmzp(ngrid)
           ind = ind + 1
           s1 = scratch%scr2(ind)
           s2 = scratch%scr1(ind)
           s3 = scratch%vt3dh(ind)
           scratch%scr2(ind) = turb_g(ngrid)%hkm(k,i,j)
           scratch%scr1(ind) = turb_g(ngrid)%vkm(k,i,j)
           scratch%vt3dh(ind) = turb_g(ngrid)%vkh(k,i,j)
           !! also for vt3di = K(tke) ?????    22 March 02
           !!         scratch%vt3di(ind) = turb_g(ngrid)%vke(k,i,j)
           turb_g(ngrid)%hkm(k,i,j) = s1
           turb_g(ngrid)%vkm(k,i,j) = s2
           turb_g(ngrid)%vkh(k,i,j) = s3
        enddo
     enddo
  enddo

  if (if_adap == 0) then

     call diffvel(mzp,mxp,myp,ia,iz,ja,jz,jdim,ia_1,ja_1             &
          ,ia1,ja1,iz_1,jz_1,iz1,jz1,izu,jzv,idiffk(ngrid)             &
          ,basic_g(ngrid)%up    (:,:,:) ,basic_g(ngrid)%vp    (:,:,:)  &
          ,basic_g(ngrid)%wp    (:,:,:) ,tend%ut              (:)      &
          ,tend%vt              (:)     ,tend%wt              (:)      &
          ,scratch%vt3da        (:)     ,scratch%vt3db        (:)      &
          ,scratch%vt3dc        (:)     ,scratch%vt3dd        (:)      &
          ,scratch%vt3de        (:)     ,scratch%vt3df        (:)      &
          ,scratch%vt3dg        (:)     ,scratch%vt3dj        (:)      &
          ,scratch%vt3dk        (:)     ,scratch%vt3dl        (:)      &
          ,scratch%vt3dm        (:)     ,scratch%vt3dn        (:)      &
          ,scratch%vt3do        (:)     ,grid_g(ngrid)%rtgu   (:,:)    &
          ,grid_g(ngrid)%rtgv   (:,:)   ,grid_g(ngrid)%rtgt   (:,:)    &
          ,turb_g(ngrid)%sflux_u(:,:)   ,turb_g(ngrid)%sflux_v(:,:)    &
          ,turb_g(ngrid)%sflux_w(:,:)   ,basic_g(ngrid)%dn0   (:,:,:)  &
          ,basic_g(ngrid)%dn0u  (:,:,:) ,basic_g(ngrid)%dn0v  (:,:,:)  &
          ,scratch%scr1         (:)     ,scratch%scr2         (:),ibcon,mynum)

  else

     call diffvel_adap(mzp,mxp,myp,ia,iz,ja,jz,jdim                  &
          ,iz1,jz1,izu,jzv,idiffk(ngrid)                               &
          ,basic_g(ngrid)%up    (:,:,:) ,basic_g(ngrid)%vp    (:,:,:)  &
          ,basic_g(ngrid)%wp    (:,:,:) ,tend%ut              (:)      &
          ,tend%vt              (:)     ,tend%wt              (:)      &
          ,scratch%vt3da        (:)     ,scratch%vt3db        (:)      &
          ,scratch%vt3dc        (:)     ,scratch%vt3dd        (:)      &
          ,scratch%vt3de        (:)     ,scratch%vt3df        (:)      &
          ,scratch%vt3dg        (:)     ,scratch%vt3dj        (:)      &
          ,scratch%vt3dk        (:)     ,scratch%vt3dl        (:)      &
          ,scratch%vt3dm        (:)     ,scratch%vt3dn        (:)      &
          ,scratch%vt3do        (:)     ,grid_g(ngrid)%aru    (:,:,:)  &
          ,grid_g(ngrid)%arv    (:,:,:) ,grid_g(ngrid)%arw    (:,:,:)  &
          ,grid_g(ngrid)%volu   (:,:,:) ,grid_g(ngrid)%volv   (:,:,:)  &
          ,grid_g(ngrid)%volw   (:,:,:) ,grid_g(ngrid)%lpu    (:,:)    &
          ,grid_g(ngrid)%lpv    (:,:)   ,grid_g(ngrid)%lpw    (:,:)    &
          ,turb_g(ngrid)%sflux_u(:,:)   ,turb_g(ngrid)%sflux_v(:,:)    &
          ,turb_g(ngrid)%sflux_w(:,:)   ,basic_g(ngrid)%dn0   (:,:,:)  &
          ,basic_g(ngrid)%dn0u  (:,:,:) ,basic_g(ngrid)%dn0v  (:,:,:)  &
          ,scratch%scr1         (:)     ,scratch%scr2         (:)      &
          ,grid_g(ngrid)%topma  (:,:)   ,ibcon,mynum)

  endif


  ! Convert momentum K's to scalar K's, if necessary

!--(DMK-CCATT)---------------------------------------------------------
  !-ml/srf - for new turb scheme
  if (idiffk(ngrid) <= 3 .or. idiffk(ngrid) == 7) then
!--(DMK-original)------------------------------------------------------
!  if (idiffk(ngrid) <= 3) then
!--(DMK-CCATT-END)-----------------------------------------------------

     do ind = 1,mxyzp
        scratch%scr2(ind) = scratch%scr2(ind) * xkhkm(ngrid)
     enddo
  elseif (idiffk(ngrid) == 4) then
     do ind = 1,mxyzp
        scratch%vt3di(ind) = 2. * scratch%scr1(ind)
     enddo
  endif


  !- CATT
!--(DMK-CCATT)---------------------------------------------------------
  if (CATT == 1 .and. CHEMISTRY >= 0) then
     ind = 0
     do j = 1,mmyp(ngrid)
        do i = 1,mmxp(ngrid)
           do k = 1,mmzp(ngrid)
              ind = ind + 1

              s1 = scratch%scr3(ind)
              scratch%scr3(ind)  = turb_s(ngrid)%hksc(k,i,j)
              turb_s(ngrid)%hksc(k,i,j) = s1

              ! salva o atual coef para o proximo passo no tempo.
           enddo
        enddo
     enddo


!--(DMK-CCATT)---------------------------------------------------------
     if (idiffk(ngrid) .le. 3 .or. idiffk(ngrid) == 7) then
!--(DMK-original)------------------------------------------------------
!     if (idiffk(ngrid) .le. 3) then
!--(DMK-CCATT-END)-----------------------------------------------------

        ind = 0
        do j = 1,mmyp(ngrid)
           do i = 1,mmxp(ngrid)
              do k = 1,mmzp(ngrid)
                 ind = ind + 1
                 scratch%scr3(ind)  =  scratch%scr3(ind)  * xkhkm(ngrid)
              enddo
           enddo
        enddo

     endif
  endif

  do n = 1,num_scalar(ngrid)

     scalarp => scalar_tab(n,ngrid)%var_p
     scalart => scalar_tab(n,ngrid)%var_t

     call azero(mxp*myp,scratch%vt2da(:))

     if (nstbot == 1) then
        if (scalar_tab(n,ngrid)%name == 'THP') then
           call atob(mxp*myp,turb_g(ngrid)%sflux_t(:,:),scratch%vt2da(:))

!--(DMK-CCATT)---------------------------------------------------------
           ! CATT
           if (CATT == 1 .and. nnqparm(ngrid) == 2) then
              !---------------------------------------
              !srf- Large Scale Forcing for GRELL CUPAR
              lsfcupar_p => grell_g_sh(ngrid)%lsfth(1,1,1)
              !---------------------------------------
           endif
!--(DMK-original)------------------------------------------------------
!           ! Large Scale Forcing for GRELL CUPAR not used
! !!$           ! CATT
! !!$           if (CATT == 1 .and. nnqparm(ngrid) == 2) then
! !!$              !---------------------------------------
! !!$              !srf- Large Scale Forcing for GRELL CUPAR
! !!$              lsfcupar_p => grell_g_sh(ngrid)%lsfth(1,1,1)
! !!$              !---------------------------------------
! !!$           endif
!--(DMK-CCATT-END)-----------------------------------------------------

        elseif (scalar_tab(n,ngrid)%name == 'RTP') then
           call atob(mxp*myp,turb_g(ngrid)%sflux_r(:,:),scratch%vt2da(:))

!--(DMK-CCATT)---------------------------------------------------------
           ! CATT
           if (CATT == 1 .and. nnqparm(ngrid) == 2) then
              !---------------------------------------
              !srf- Large Scale Forcing for GRELL CUPAR
              lsfcupar_p => grell_g_sh(ngrid)%lsfrt(1,1,1)
              !---------------------------------------
           endif
!--(DMK-original)------------------------------------------------------
!           ! Large Scale Forcing for GRELL CUPAR not used
! !!$           ! CATT
! !!$           if (CATT == 1 .and. nnqparm(ngrid) == 2) then
! !!$              !---------------------------------------
! !!$              !srf- Large Scale Forcing for GRELL CUPAR
! !!$              lsfcupar_p => grell_g_sh(ngrid)%lsfrt(1,1,1)
! !!$              !---------------------------------------
! !!$           endif
!--(DMK-CCATT-END)-----------------------------------------------------

        endif
     endif

     ! 3/10/01 - Define ksf below, the "K scalar flag", to let subroutine diffsclr
     ! know which vertical K is being passed to it.  If diffsclr sees that it's
     ! a different K from the previous one, diffsclr will re-compute the tridiff
     ! matrix coefficients.  In order to use vertical scalar K's other than
     ! vt3dh and vt3di, use ksf = 3, ksf = 4, etc. for each different K.

     !_STC..................................................
     !_STC Corrections to account for the new idiffk options
     !_STC for E-l and E-eps closure. Isotropy hypothesis.
     !_STC (S. Trini Castelli)
     !_STC..................................................

     if (scalar_tab(n,ngrid)%name == 'TKEP') then
        vkh_p => scratch%vt3di(1)
        hkh_p => scratch%scr2(1)

!--(DMK-CCATT)---------------------------------------------------------
        !-ml/srf - for new turb scheme        
        if (idiffk(ngrid) >= 4 .and. idiffk(ngrid) /= 7) hkh_p => scratch%vt3di(1)
!--(DMK-original)------------------------------------------------------     
!        if (idiffk(ngrid) >= 4) hkh_p => scratch%vt3di(1)
!--(DMK-CCATT-END)-----------------------------------------------------

        ksf = 1
     elseif (scalar_tab(n,ngrid)%name == 'EPSP') then
        vkh_p => scratch%vt3di(1)
        hkh_p => scratch%scr2(1)

!--(DMK-CCATT)---------------------------------------------------------
        !-ml/srf - for new turb scheme   
        if (idiffk(ngrid) >= 4 .and. idiffk(ngrid) /= 7) hkh_p => scratch%vt3di(1)
!--(DMK-original)------------------------------------------------------     
!        if (idiffk(ngrid) >= 4)  hkh_p => scratch%vt3di(1)
!--(DMK-CCATT-END)-----------------------------------------------------

        ksf = 3
        ! Convert Ktke to Keps; it will be converted back after use below
        call ae1t0 (mxyzp,vkh_p,vkh_p,ALF_EPS/ALF_TKE)
        call ae1t0 (mxyzp,hkh_p,hkh_p,ALF_EPS/ALF_TKE)
     else
        vkh_p => scratch%vt3dh(1)
        hkh_p => scratch%scr2(1)

!--(DMK-CCATT)---------------------------------------------------------
        !-ml/srf - for new turb scheme        
        if (idiffk(ngrid) >= 4 .and. idiffk(ngrid) /= 7) hkh_p => scratch%vt3dh(1)
!--(DMK-original)------------------------------------------------------     
!        if (idiffk(ngrid) >= 4) hkh_p => scratch%vt3dh(1)
!--(DMK-CCATT-END)-----------------------------------------------------

        ksf = 2
     endif

!--(DMK-CCATT)---------------------------------------------------------
     ! CATT
     if (CATT == 1 .and. CHEMISTRY >= 0) then
        if(n > (num_scalar(ngrid) - (NADDSC + NSPECIES_TRANSPORTED))) then
           if (idiffk(ngrid) < 4 .or.  idiffk(ngrid) == 7) then
              hkh_p => scratch%scr3(1)  !ivt3dp
           endif
        endif
     endif
!--(DMK-original)------------------------------------------------------     
!     if (CATT == 1) then
!        !srf----------------- Hor. Diffusion Coef for tracers
!        if(n > (num_scalar(ngrid) - NADDSC)) then ! if(n .gt. nscl-NADDSC)
!           if (idiffk(ngrid) /= 4) then
!              hkh_p => scratch%scr3(1)  !hkh_pscr3 => scratch%vt3dp(1) !ivt3dp
!           endif
!        endif
!        !srf--------------------
!     endif
!--(DMK-CCATT-END)-----------------------------------------------------

     if (if_adap == 0) then

        call diffsclr(mzp,mxp,myp,ia,iz,ja,jz,jdim                  &
             ,ia_1,ja_1,ia1,ja1,iz_1,jz_1,iz1,jz1,n,ksf               &
             ,scalarp,scalart            ,scratch%vt3da(:)            &
             ,scratch%vt3db      (:)     ,scratch%vt3df      (:)      &
             ,scratch%vt3dg      (:)     ,scratch%vt3dj      (:)      &
             ,scratch%vt3dk      (:)     ,scratch%vt3do      (:)      &
             ,scratch%vt3dc      (:)     ,scratch%vt3dd      (:)      &
             ,scratch%vt3dl      (:)     ,scratch%vt3dm      (:)      &
             ,scratch%vt2db      (:)     ,grid_g(ngrid)%rtgt (:,:)    &
             ,scratch%vt2da      (:)     ,basic_g(ngrid)%dn0 (:,:,:)  &
             ,vkh_p                      ,hkh_p                       &!)
             !srf - large and subgrid scale for GRELL CUPAR (CATT)
!--(DMK-CCATT)---------------------------------------------------------
             ,lsfcupar_p,n)
!--(DMK-original)------------------------------------------------------     
!        ! large and subgrid scale for GRELL CUPAR (CATT) not used
!        !       ,lsfcupar_p,n)
!--(DMK-CCATT-END)-----------------------------------------------------
     else

        call diffsclr_adap(mzp,mxp,myp,ia,iz,ja,jz,jdim,n,ksf      &
             ,grid_g(ngrid)%lpw(:,:)     ,scalarp                    &
             ,scalart                    ,scratch%vt3da     (:)      &
             ,scratch%vt3dc      (:)     ,scratch%vt3df     (:)      &
             ,scratch%vt3dg      (:)     ,scratch%vt3dj     (:)      &
             ,scratch%vt3dk      (:)     ,scratch%vt3dl     (:)      &
             ,scratch%vt3dm      (:)     ,scratch%vt3do     (:)      &
             ,scratch%vt2da      (:)     ,scratch%vt2db     (:)      &
             ,basic_g(ngrid)%dn0 (:,:,:) ,vkh_p                      &
             ,hkh_p                      ,grid_g(ngrid)%aru (:,:,:)  &
             ,grid_g(ngrid)%arv  (:,:,:) ,grid_g(ngrid)%arw (:,:,:)  &
             ,grid_g(ngrid)%volt (:,:,:) ,scratch%vt3db     (:)      &
             ,grid_g(ngrid)%dxu  (:,:)   ,grid_g(ngrid)%dyv (:,:)    &
             ,grid_g(ngrid)%topma(:,:)                               &
             !srf - large and subgrid scale for GRELL CUPAR (CATT)
!--(DMK-CCATT)---------------------------------------------------------
             ,lsfcupar_p,n)
!--(DMK-original)------------------------------------------------------     
!        ! CATT not ready for Shaved-ETA
!        !     ,lsfcupar_p,n)
!--(DMK-CCATT-END)-----------------------------------------------------
     endif

     if (scalar_tab(n,ngrid)%name == 'EPSP') then
        call ae1t0 (mxyzp,vkh_p,vkh_p,ALF_TKE/ALF_EPS)
        call ae1t0 (mxyzp,hkh_p,hkh_p,ALF_TKE/ALF_EPS)
     endif

  enddo

  return
end subroutine diffuse

!     ******************************************************************

subroutine strain(m1,m2,m3,ia,iz,ja,jz,ia_1,ja_1,iz1,jz1  &
     ,jd,up,vp,wp,vt3da,vt3db,vt3dc,vt3dd,vt3de  &
     ,vt3df,vt3dg,vt3dh,vt3di,vt3dn,scr2,idiffk)

  implicit none

  integer, intent(in)  :: m1   &
                        , m2   &
                        , m3   &
                        , ia   &
                        , iz   &
                        , ja   &
                        , jz   &
                        , ia_1 &
                        , ja_1 &
                        , iz1  &
                        , jz1  &
                        , jd   &
                        , idiffk

  real, dimension(m1,m2,m3), intent(in) :: up,vp,wp

  real, dimension(m1,m2,m3), intent(inout):: vt3da,      &
                                             vt3db,      &
                                             vt3dc,      &
                                             vt3dd,      &
                                             vt3de,      &
                                             vt3df,      &
                                             vt3dg,      &
                                             vt3dh,      &
                                             vt3di,      &
                                             vt3dn,      &
                                             scr2

  !local variables:
  integer              :: i,j,k


  call grad(m1, m2, m3, ia  , iz1, ja  , jz , up, vt3da, 'XDIR', 'UPNT')
  call grad(m1, m2, m3, ia_1, iz , ja_1, jz , vp, vt3db, 'XDIR', 'VPNT')
  call grad(m1, m2, m3, ia_1, iz , ja  , jz , wp, vt3df, 'XDIR', 'WPNT')

  call grad(m1, m2, m3, ia_1, iz , ja_1, jz , up, vt3dn, 'YDIR', 'UPNT')
  call grad(m1, m2, m3, ia  , iz , ja  , jz1, vp, vt3dc, 'YDIR', 'VPNT')
  call grad(m1, m2, m3, ia  , iz , ja_1, jz , wp, vt3dg, 'YDIR', 'WPNT')

  call grad(m1, m2, m3, ia_1, iz , ja  , jz , up, vt3dd, 'ZDIR', 'UPNT')
  call grad(m1, m2, m3, ia  , iz , ja_1, jz , vp, vt3de, 'ZDIR', 'VPNT')

!--(DMK-CCATT)---------------------------------------------------------
  !ml/srf- for new turb scheme
  if(idiffk >= 3 .and. idiffk /= 7)then 
!--(DMK-original)------------------------------------------------------
!  if(idiffk.ge.3)then
!--(DMK-CCATT-END)-----------------------------------------------------

     call grad(m1,m2,m3,ia,iz,ja,jz,wp,scr2,'ZDIR','WPNT')
  endif

!--(DMK-CCATT)---------------------------------------------------------
  !ml/srf- for new turb scheme
  if (idiffk <= 2 .or. idiffk == 7) then
!--(DMK-original)------------------------------------------------------
!  if (idiffk .le. 2) then
!--(DMK-CCATT-END)-----------------------------------------------------

     do j = ja,jz
        do i = ia,iz
           do k = 2,m1-1
              vt3dh(k,i,j) =2. * (vt3da(k,i,j) * vt3da(k,i,j)  &
                   + vt3dc(k,i,j) * vt3dc(k,i,j))  &
                   + .0625 * (vt3db(k,i,j) + vt3db(k,i-1,j)  &
                   + vt3db(k,i,j-jd) + vt3db(k,i-1,j-jd)  &
                   + vt3dn(k,i,j) + vt3dn(k,i-1,j)  &
                   + vt3dn(k,i,j-jd) + vt3dn(k,i-1,j-jd)) ** 2
              vt3di(k,i,j) = .0625 * ((vt3dd(k,i,j) + vt3dd(k-1,i,j)  &
                   + vt3dd(k,i-1,j) + vt3dd(k-1,i-1,j)) ** 2  &
                   + (vt3de(k,i,j) + vt3de(k-1,i,j)  &
                   + vt3de(k,i,j-jd) + vt3de(k-1,i,j-jd)) ** 2)
           enddo
        enddo
     enddo
  else
     do j = ja,jz
        do i = ia,iz
           do k = 2,m1-1
              vt3da(k,i,j) = 2. * vt3da(k,i,j)
              vt3dc(k,i,j) = 2. * vt3dc(k,i,j)
              scr2(k,i,j) = 2. * scr2(k,i,j)
              vt3db(k,i,j) = vt3db(k,i,j) + vt3dn(k,i,j)
              vt3dn(k,i,j) = vt3db(k,i,j)
              vt3dd(k,i,j) = vt3dd(k,i,j) + vt3df(k,i,j)
              vt3de(k,i,j) = vt3de(k,i,j) + vt3dg(k,i,j)
              vt3di(k,i,j) = 0.333333  &
                   * (vt3da(k,i,j) + vt3dc(k,i,j) + scr2(k,i,j))
           enddo
        enddo

        do k = 2,m1-1
           vt3da(k,iz1,j) = 2. * vt3da(k,iz1,j)
           vt3db(k,ia_1,j) = vt3db(k,ia_1,j) + vt3dn(k,ia_1,j)
           vt3dn(k,ia_1,j) = vt3db(k,ia_1,j)
           vt3dd(k,ia_1,j) = vt3dd(k,ia_1,j) + vt3df(k,ia_1,j)
        enddo
     enddo

     do i = ia_1,iz
        do k = 2,m1-1
           vt3dc(k,i,jz1) = 2. * vt3dc(k,i,jz1)
           vt3db(k,i,ja_1) = vt3db(k,i,ja_1) + vt3dn(k,i,ja_1)
           vt3dn(k,i,ja_1) = vt3db(k,i,ja_1)
           vt3de(k,i,ja_1) = vt3de(k,i,ja_1) + vt3dg(k,i,ja_1)
        enddo
     enddo

     do j = ja,jz
        do i = ia,iz
           do k = 2,m1-1
              vt3dh(k,i,j) = .5 * (  &
                   (vt3da(k,i,j) - vt3di(k,i,j)) ** 2  &
                   + (vt3dc(k,i,j) - vt3di(k,i,j)) ** 2  &
                   + ( scr2(k,i,j) - vt3di(k,i,j)) ** 2)  &
                   + .0625 * ((vt3db(k,i,j) + vt3db(k,i-1,j)  &
                   + vt3db(k,i,j-jd) + vt3db(k,i-1,j-jd)) ** 2  &
                   + (vt3dd(k,i,j) + vt3dd(k,i-1,j)  &
                   + vt3dd(k-1,i,j) + vt3dd(k-1,i-1,j)) ** 2  &
                   + (vt3de(k,i,j) + vt3de(k-1,i,j)  &
                   + vt3de(k,i,j-jd) + vt3de(k-1,i,j-jd)) ** 2)
              vt3di(k,i,j) = vt3dh(k,i,j)
           enddo
        enddo
     enddo
  endif

  return
end subroutine strain

!     *****************************************************************

subroutine bruvais(m1,m2,m3,ia,iz,ja,jz,theta,rtp,rv,rcp,pp,pi0,en2  &
     ,rtgt,lpw)

  use mem_scratch, only : vctr11,    &     !INTENT(INOUT)
                          vctr12,    &     !INTENT(INOUT)
                          vctr32,    &     !INTENT(INOUT)
                          vctr1,     &     !INTENT(INOUT)
                          vctr2,     &     !INTENT(INOUT)
                          vctr3,     &     !INTENT(INOUT)
                          vctr4            !INTENT(INOUT)

  use micphys, only     : level,     &     !INTENT(in)
                          ipris,     &     !INTENT(in)
                          isnow,     &     !INTENT(in)
                          igraup,    &     !INTENT(in)
                          iaggr,     &     !INTENT(in)
                          ihail            !INTENT(in)

  use mem_grid, only    : zt,        &     !INTENT(in)
                          nzp,       &     !INTENT(in)
                          nz,        &     !INTENT(in)
                          nzpmax           !INTENT(in)

  use rconstants, only  : alvl,      &     !INTENT(in)
                          rgas,      &     !INTENT(in)
                          ep,        &     !INTENT(in)
                          cp,        &     !INTENT(in)
                          alvi,      &     !INTENT(in)
                          p00,       &     !INTENT(in)
                          cpor,      &     !INTENT(in)
                          g                !INTENT(in)

  implicit none

  integer, INTENT(IN) :: m1  &
                       , m2  &
                       , m3  &
                       , ia  &
                       , iz  &
                       , ja  &
                       , jz

  real, dimension(m1,m2,m3), INTENT(in)    :: theta,   &
                                              pp,      &
                                              pi0,     &
                                              rtp,     &
                                              rcp,     &
                                              rv


  real, dimension(m1,m2,m3), INTENT(INOUT)   :: en2

  real, dimension(m2,m3), INTENT(IN)         :: rtgt

  integer, dimension(m2,m3), INTENT(IN)      :: lpw

  !local variables
  integer :: i,j,k,iweten,iwdiffk,ki,k2,k1

  real :: c1,c2,c3,ci1,ci2,ci3,rvlsi,rvii
  real, dimension(nzpmax) :: pi,temp,prt,rvls,rc
  ! **(JP)** fatora expressoes logicas para fora dos lacos
  logical :: log1, log2, log3, log4


  !     calculate brunt-vaisalla frequency squared (en2)

  iweten = 1

  iwdiffk = 0
  c1 = alvl / rgas
  c2 = ep * alvl ** 2 / (cp * rgas)
  c3 = alvl / cp
  ci1 = alvi / rgas
  ci2 = ep * alvi ** 2 / (cp * rgas)
  ci3 = alvi / cp

  ! **(JP)** fatora expressoes logicas para fora dos lacos
  log1 = level .ge. 1
  log2 = level .ge. 2 .and. iweten .eq. 1
  log3 = (ipris  .ge. 1 .or. isnow .ge. 1 .or.  &
       igraup .ge. 1 .or. iaggr .ge. 1 .or.  &
       ihail  .ge. 1) .and. level .eq. 3
  log4 = level .eq. 3
  ! **(JP)** fim modificacao

  !     calculate potential temperature profile

  do j = ja,jz
     do i = ia,iz

        k2=lpw(i,j)
        k1=k2-1

        do k = k1,m1
           vctr11(k) = theta(k,i,j)
           vctr12(k) = theta(k,i,j)
           vctr32(k) = 0.
        enddo
        ! **(JP)** fatora expressoes logicas para fora dos lacos
        if (log1) then
        ! **(JP)** fim modificacao
           do k = k1,m1
              vctr12(k) = vctr11(k) * (1. + .61 * rv(k,i,j))
              vctr32(k) = (rtp(k,i,j) - rv(k,i,j))
           enddo
        endif

        !     check for saturation if level is 2 or greater.
        ! **(JP)** fatora expressoes logicas para fora dos lacos
        if (log2) then
        ! **(JP)** fim modificacao
           do k = k1,m1
              pi(k) = (pp(k,i,j) + pi0(k,i,j)) / cp
              temp(k) = theta(k,i,j) * pi(k)
              prt(k) = p00 * pi(k) ** cpor
           enddo
           call mrsl(m1,prt(:),temp(:),rvls(:))
           do k = k2-1,m1
              vctr2(k) = c1
              vctr3(k) = c2
              vctr4(k) = c3
           enddo
           ki = m1 + 1

           !     if any ice phase microphysics are activated ....
           ! **(JP)** fatora expressoes logicas para fora dos lacos
           if (log3) then
           ! **(JP)** fim modificacao
              !  find level of -20 c.  assume ice saturation above this
              !   level.

              do k = k1,m1
                 if (temp(k) .le. 253.16) then
                    ki = k
                    go to 10
                 endif
              enddo
              ki = m1 + 1
10            continue
              call mrsi(m1-ki+1,prt(ki),temp(ki),rvls(ki))
!-srf 19/03/2005
!bug: Subscript out of range for array prt
!    subscript=0, lower bound=1, upper bound=132, dimension=1
! tempc < 253 sobre parte da Antartica.
!              ki = max(ki,k2)                                    !solucao 1
!              call mrsi(1,prt(ki-1),temp(ki-1),rvlsi)
               if(ki > 1) call mrsi(1,prt(ki-1),temp(ki-1),rvlsi) ! solucao 2
!srf


              do k = ki,m1
                 vctr2(k) = ci1
                 vctr3(k) = ci2
                 vctr4(k) = ci3
              enddo
           endif

           ! **(JP)** fatora expressoes logicas para fora dos lacos
           if (log4) then
           ! **(JP)** fim modificacao
              do k = k1,m1
                 rc(k) = rcp(k,i,j)
              enddo
           else
              do k = k1,m1
                 rc(k) = max(rv(k,i,j) / rvls(k) - .999,0.)
              enddo
           endif

        endif

        do k = k2,m1-1
           vctr1(k) = g / ((zt(k+1) - zt(k-1)) * rtgt(i,j))
        enddo

        ! **(JP)** fatora expressoes logicas para fora dos lacos
        if (log2) then
        ! **(JP)** fim modificacao
           do k = k2,m1-1
              if (rc(k) .gt. 0.) then
                 rvii = rvls(k-1)
                 if (k .eq. ki) rvii = rvlsi
                 en2(k,i,j) = vctr1(k) * (  &
                      (1. + vctr2(k) * rvls(k) / temp(k))  &
                      / (1. + vctr3(k) * rvls(k) / temp(k) ** 2)  &
                      * ((vctr11(k+1) - vctr11(k-1)) / vctr11(k)  &
                      + vctr4(k) / temp(k) * (rvls(k+1) - rvii))  &
                      - (rtp(k+1,i,j) - rtp(k-1,i,j)))
              else
                 en2(k,i,j) = vctr1(k)*((vctr12(k+1)-vctr12(k-1))  &
                      / vctr12(k) - (vctr32(k+1) - vctr32(k-1)))
              endif

           enddo
        else
           do k = k2,m1-1
              en2(k,i,j) = vctr1(k) * ((vctr12(k+1)-vctr12(k-1))  &
                   / vctr12(k) - (vctr32(k+1) - vctr32(k-1)))
           enddo
        endif
        ! **(JP)** remove para fora do laco, permitindo vetorizacao
        !en2(k1,i,j) = en2(k2,i,j)
        !en2(nzp,i,j)=en2(nz,i,j)
        ! **(JP)** fim da modificacao

     enddo
  enddo

  ! **(JP)** removido de dentro do laco
  do j = ja,jz
     do i = ia,iz
        en2(lpw(i,j)-1,i,j)=en2(lpw(i,j),i,j)
     end do
  end do
  do j = ja,jz
     do i = ia,iz
        en2(nzp,i,j)=en2(nz,i,j)
     end do
  end do
  ! **(JP)** fim da modificacao

  return
end subroutine bruvais

!     *****************************************************************

subroutine mxdefm(m1,m2,m3,ia,iz,ja,jz,ibcon,jd  &
     ,vt3dh,vt3di,vt3dj,vt3dk,scr1,scr2,dn0,rtgt,dxt,dyt,lpw,mynum)

  !     +-------------------------------------------------------------+
  !     \   this routine calculates the mixing coefficients with a    \
  !     \     smagorinsky-type deformational based k with an optional \
  !     \     unstable brunt-vaisala enhancement and an optional      \
  !     \     richardson number modification.                         \
  !     +-------------------------------------------------------------+

  use mem_scratch , only   : vctr1,    &     !INTENT(INOUT)
                             vctr2           !INTENT(INOUT)

  use mem_grid, only       : ngrid,    &     !INTENT(in)
                             zm,       &     !INTENT(in)
                             zt              !INTENT(in)

  use mem_turb, only       : csx,      &     !INTENT(in)
                             idiffk,   &     !INTENT(in)
                             rmin,     &     !INTENT(out)
                             rmax,     &     !INTENT(out)
                             zkhkm,    &     !INTENT(in)
                             csz,      &     !INTENT(in)
                             akmin           !INTENT(in)

  use rconstants, only     : vonk
 
  !-srf para diferente difusao numerica >
  use extras, only:   extra2d,NA_EXTRA2D ! intent(in)
  use catt_start, only: CATT           ! intent(in)

  implicit none

  integer, intent(in) :: m1   &     !INTENT(in)
                       , m2   &     !INTENT(in)
                       , m3   &     !INTENT(in)
                       , ia   &     !INTENT(in)
                       , iz   &     !INTENT(in)
                       , ja   &     !INTENT(in)
                       , jz         !INTENT(in)

  integer, intent(in) :: ibcon, jd, mynum   !- EHE -> nao sao usadas!!!

  real, dimension(m1,m2,m3), INTENT(INOUT) :: vt3dh    &
                                            , vt3dk    &
                                            , scr1     &
                                            , scr2

  real, dimension(m1,m2,m3), INTENT(IN)    :: dn0      &
                                            , vt3di    &
                                            , vt3dj

  real, dimension(m2,m3), INTENT(IN) :: rtgt,dxt

  real, dimension(m2,m3), INTENT(IN) :: dyt  !- EHE -> nao e' usada!!!

  integer, dimension(m2,m3), INTENT(IN) :: lpw


  !local variables:
  integer :: i,j,k,irich,ienfl

  real :: csx2,sq300,enfl,rchmax,c1,c2,c3,c4,akm,ambda,vkz2

  irich = 1
  ienfl = 1

  csx2 = csx(ngrid) * csx(ngrid)
  sq300 = 90000.
  if (idiffk(ngrid) .eq. 2 .or. idiffk(ngrid) .eq. 3) then
     rmin = -100.
     rmax = 1. / zkhkm(ngrid)
     do j = ja,jz
        do i = ia,iz
           do k = lpw(i,j),m1-1
              vt3dk(k,i,j) = max(min(vt3dj(k,i,j)  &
                   / max(vt3di(k,i,j),1.e-15),rmax),rmin)
           enddo
        enddo
     enddo
     enfl = float(ienfl)
     rchmax = 1.0 + 9.0 * float(irich)

!--(DMK-CCATT)---------------------------------------------------------
     !mudanca rams60 nova versao
     do k = 2, m1
!--(DMK-original)------------------------------------------------------
!     do k = lpw(i,j),m1
!--(DMK-CCATT)---------------------------------------------------------

        vctr1(k) = csz(ngrid) * (zm(k) - zm(k-1))
        vctr2(k) = vctr1(k) * vctr1(k)
     enddo
  endif

!--(DMK-CCATT)---------------------------------------------------------
  !ml/srf - for new turb scheme
  if (idiffk(ngrid) == 1 .or. idiffk(ngrid) == 7) then

    if(CATT==1 .and. akmin(ngrid) < 0.) then 
      !-srf: testa numero de extras:
         if(NA_EXTRA2D < 5) stop 'NAEXTRA2d must be at least 5' 
      !srf
   	do j = ja,jz
   	    do i = ia,iz
   	       c2 = 1.0 / (dxt(i,j) * dxt(i,j))
   	       c3 = csx2 * c2
   	       akm = abs(akmin(ngrid)) * 0.075 * c2 ** (0.666667)
   	       !----
   	       !srf-define diferentes AKMINs para melhorar estabilidade 
   	       !srf-sobre os Andes
   	       akm=extra2d(5,ngrid)%d2(i,j)*akm
   	       !----
   	       do k = lpw(i,j),m1-1
   		  scr2(k,i,j) = dn0(k,i,j)  &
   		       * max(akm,c3*sqrt(vt3dh(k,i,j)))
   	       enddo
   	    enddo
   	 enddo
    else
   	do j = ja,jz
   	    do i = ia,iz
   	       c2 = 1.0 / (dxt(i,j) * dxt(i,j))
   	       c3 = csx2 * c2
   	       akm = akmin(ngrid) * 0.075 * c2 ** (0.666667)
   	       do k = lpw(i,j),m1-1
   		  scr2(k,i,j) = dn0(k,i,j)  &
   		       * max(akm,c3*sqrt(vt3dh(k,i,j)))
   	       enddo
   	    enddo
   	 enddo
    endif
    
  elseif (idiffk(ngrid) .eq. 2) then
     do j = ja,jz
        do i = ia,iz
           c1 = rtgt(i,j) * rtgt(i,j)
           c2 = 1.0 / (dxt(i,j) * dxt(i,j))
           c3 = csx2 * c2
           akm = akmin(ngrid) * 0.075 * c2 ** (0.666667)
           c4 = vonk * vonk * c1
           do k = lpw(i,j),m1-1
              ! old csz*dz len  scr1(k,i,j) = dn0(k,i,j) * c1 * vctr2(k)

              ! asymptotic vertical scale length from bjorn with modifications:
              ! c3 is (csx * dx)^2, c1*vctr2(k) is (csz * dz)^2, sq300 is the square
              ! of 300 meters (used as a limit for horizontal grid spacing influence
              ! on vertical scale length), ambda is (asymptotic_vertical_length_scale)^2,
              ! and vkz2 is (vonk * height_above_surface)^2.

              ambda = max(c1 * vctr2(k),min(sq300,c3))
              vkz2 = c4 * zt(k) * zt(k)
              scr1(k,i,j) = dn0(k,i,j) * vkz2 / (vkz2 / ambda + 1)  &

              * (sqrt(vt3di(k,i,j))  &
                   + enfl * sqrt(max(0.,-vt3dj(k,i,j))))*min(rchmax  &
                   ,sqrt(max(0.,(1.-zkhkm(ngrid)*vt3dk(k,i,j)))))

              scr2(k,i,j) = dn0(k,i,j)  &
                   * max(akm,c3*sqrt(vt3dh(k,i,j)))
              vt3dh(k,i,j) = scr1(k,i,j) * zkhkm(ngrid)

           enddo
        enddo
     enddo
     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     !     call friclyr(nzp,nxp,nyp,a(iscr1),a(iustarl),a(itstarl)
     !    +    ,a(iustarw),a(itstarw),a(ipctlnd),a(itheta),a(irtgt))
     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  elseif (idiffk(ngrid) .eq. 3) then
     do j = ja,jz
        do i = ia,iz
           c1 = rtgt(i,j) * rtgt(i,j)
           do k = lpw(i,j),m1-1
              scr1(k,i,j) = dn0(k,i,j) * c1 * vctr2(k)  &
                   * (sqrt(vt3dh(k,i,j))  &
                   + enfl * sqrt(max(0.,-vt3dj(k,i,j))))*min(rchmax  &
                   ,sqrt(max(0.,(1.-zkhkm(ngrid)*vt3dk(k,i,j)))))
              scr2(k,i,j) = scr1(k,i,j)
              vt3dh(k,i,j) = scr1(k,i,j) * zkhkm(ngrid)
           enddo
        enddo
     enddo
     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     !     call friclyr(nzp,nxp,nyp,a(iscr1),a(iustarl),a(itstarl)
     !    +    ,a(iustarw),a(itstarw),a(ipctlnd),a(itheta),a(irtgt))
     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  endif

  return
end subroutine mxdefm

!     *****************************************************************

subroutine klbnd(m1,m2,m3,ibcon,jd,akay,dn0,lpw)

  implicit none

  integer, INTENT(IN)       :: m1     &
                             , m2     &
                             , m3     &
                             , ibcon  &
                             , jd

  real, dimension(m1,m2,m3), INTENT(INOUT) :: akay

  real, dimension(m1,m2,m3), INTENT(IN)    :: dn0

  integer, dimension(m2,m3), INTENT(IN)    :: lpw

  !local variables:
  integer ::i,j,k,k2

  !     boundary conditions on a mixing coefficient

  do j = 1,m3
     do i = 1,m2
        k2=lpw(i,j)
        do k=1,lpw(i,j)-1
           akay(k,i,j) = akay(k2,i,j) * dn0(k,i,j) / dn0(k2,i,j)
        enddo
        akay(m1,i,j) = akay(m1-1,i,j) * dn0(m1,i,j)  &
             / dn0(m1-1,i,j)
     enddo
  enddo
  if (iand(ibcon,1) .ne. 0) then

     do j = 1,m3
        do k = 1,m1
           akay(k,1,j) = akay(k,2,j)
        enddo
     enddo
  endif

  if (iand(ibcon,2) .ne. 0) then
     do j = 1,m3
        do k = 1,m1
           akay(k,m2,j) = akay(k,m2-1,j)
        enddo
     enddo
  endif

  if (jd .eq. 1) then
     if (iand(ibcon,4) .ne. 0) then
        do i = 1,m2
           do k = 1,m1
              akay(k,i,1) = akay(k,i,2)
           enddo
        enddo
     endif

     if (iand(ibcon,8) .ne. 0) then
        do i = 1,m2
           do k = 1,m1
              akay(k,i,m3) = akay(k,i,m3-1)
           enddo
        enddo
     endif
  endif

  return
end subroutine klbnd


!     *****************************************************************
subroutine mxdefm_tracer(m1,m2,m3,ia,iz,ja,jz,ibcon,jd  &
     ,vt3dh,khtr,dn0,dxt,dyt,lpw,mynum)  !ALF

  !     +-------------------------------------------------------------+
  !     \   this routine calculates the mixing coefficients with a    \
  !     \     smagorinsky-type deformational based k                  \
  !     +-------------------------------------------------------------+
  !       khtr = coef. dif. horizontal para tracers
  !       frtr = fator de reducao do Akmin dos campos meteorologicos
  !
  use mem_grid, only:  &
       ngrid,ngrids                !INTENT(IN)
  use mem_turb, only:  &
       csx,            &    !INTENT(IN)
       akmin                !INTENT(IN)

  implicit none
  ! Arguments:
  integer, intent(in)                      :: m1, m2, m3, ia, iz, ja, jz, ibcon, jd, mynum
  ! ALF
  integer, dimension(m2,m3), intent(in)    :: lpw
  real, dimension(m2,m3), intent(in)       :: dxt, dyt    !dyt nao eh usada
  real, dimension(m1,m2,m3), intent(in)    :: vt3dh, dn0
  real, dimension(m1,m2,m3), intent(inout) ::  khtr

  ! local variables:
  integer :: i,j,k
  real :: csx2,c2,c3,akm!,frtr
  real,dimension(ngrids) :: frtr

  !frtr = 0.1
  frtr(1) = 1.
  if(ngrids> 1) frtr(2:ngrids)= 0.5
  
  csx2 = csx(ngrid) * csx(ngrid)

  do j = min(1, ja), max(jz, m3)
     do i = min(1, ia), max(iz, m2)

        c2 = 1.0 / (dxt(i,j) * dxt(i,j))
        c3 = csx2 * c2
        akm = frtr(ngrid) * abs(akmin(ngrid)) * 0.075 * c2 ** (0.666667)

        do k = lpw(i,j),m1-1

           khtr(k,i,j) = dn0(k,i,j)  &
                * max(akm,c3*sqrt(vt3dh(k,i,j)))

        enddo

     enddo
  enddo

  return
end subroutine mxdefm_tracer

!     ******************************************************************
!     *****************************************************************
!srf - OLAM
 subroutine bruvais_OLAM(m1,m2,m3,ia,iz,ja,jz,theta,rtp,rv,rcp,pp,pi0,en2  &
     ,rtgt,lpw)

  use mem_scratch, only : vctr11,    &     !INTENT(INOUT)
                          vctr12,    &     !INTENT(INOUT)
                          vctr32,    &     !INTENT(INOUT)
                          vctr1,     &     !INTENT(INOUT)
                          vctr2,     &     !INTENT(INOUT)
                          vctr3,     &     !INTENT(INOUT)
                          vctr4            !INTENT(INOUT)

  use micphys, only     : level,     &     !INTENT(in)
                          ipris,     &     !INTENT(in)
                          isnow,     &     !INTENT(in)
                          igraup,    &     !INTENT(in)
                          iaggr,     &     !INTENT(in)
                          ihail            !INTENT(in)

  use mem_grid, only    : zt,        &     !INTENT(in)
                          nzp,       &     !INTENT(in)
                          nz,        &     !INTENT(in)
                          nzpmax           !INTENT(in)

  use rconstants, only  : alvl,      &     !INTENT(in)
                          rgas,      &     !INTENT(in)
                          ep,        &     !INTENT(in)
                          cp,        &     !INTENT(in)
                          alvi,      &     !INTENT(in)
                          p00,       &     !INTENT(in)
                          cpor,      &     !INTENT(in)
                          g                !INTENT(in)

  implicit none

  integer, INTENT(IN) :: m1  &
                       , m2  &
                       , m3  &
                       , ia  &
                       , iz  &
                       , ja  &
                       , jz

  real, dimension(m1,m2,m3), INTENT(in)    :: theta,   &
                                              pp,      &
                                              pi0,     &
                                              rtp,     &
                                              rcp,     &
                                              rv


  real, dimension(m1,m2,m3), INTENT(INOUT)   :: en2

  real, dimension(m2,m3), INTENT(IN)         :: rtgt

  integer, dimension(m2,m3), INTENT(IN)      :: lpw

  !local variables
  integer :: i,j,k,iweten,iwdiffk,ki,k2,k1

  real :: c1,c2,c3,ci1,ci2,ci3,rvlsi,rvii
  real, dimension(nzpmax) :: pi,temp,prt,rvls,rc
  ! **(JP)** fatora expressoes logicas para fora dos lacos
  logical :: log1, log2, log3, log4
  
  !-------------------------------------------------
  !srf- 29/12/2008:  adaptado da versao OLAM3.0
  !-------------------------------------------------
  !     calculate brunt-vaisalla frequency squared (en2)

  !iweten = 1
  !
  !iwdiffk = 0
  !c1 = alvl / rgas
  !c2 = ep * alvl ** 2 / (cp * rgas)
  !c3 = alvl / cp
  !ci1 = alvi / rgas
  !ci2 = ep * alvi ** 2 / (cp * rgas)
  !ci3 = alvi / cp

  ! **(JP)** fatora expressoes logicas para fora dos lacos
  log1 = level .ge. 1
!  log2 = level .ge. 2 .and. iweten .eq. 1
!  log3 = (ipris  .ge. 1 .or. isnow .ge. 1 .or.  &
!       igraup .ge. 1 .or. iaggr .ge. 1 .or.  &
!       ihail  .ge. 1) .and. level .eq. 3
!  log4 = level .eq. 3
! **(JP)** fim modificacao

  !     calculate potential temperature profile

  do j = ja,jz
     do i = ia,iz

        k2=lpw(i,j)
        k1=k2-1

        do k = k1,m1
           vctr11(k) = theta(k,i,j)
           vctr12(k) = theta(k,i,j)
           !vctr32(k) = 0.
        enddo
        ! **(JP)** fatora expressoes logicas para fora dos lacos
        if (log1) then
        ! **(JP)** fim modificacao
           do k = k1,m1
              vctr12(k) = vctr11(k) * (1. + .61 * rv(k,i,j))
              !vctr32(k) = (rtp(k,i,j) - rv(k,i,j))
           enddo
        endif

        !     check for saturation if level is 2 or greater.
        ! **(JP)** fatora expressoes logicas para fora dos lacos
        !if (log2) then
        ! **(JP)** fim modificacao
        !   do k = k1,m1
        !      pi(k) = (pp(k,i,j) + pi0(k,i,j)) / cp
        !      temp(k) = theta(k,i,j) * pi(k)
        !      prt(k) = p00 * pi(k) ** cpor
        !   enddo
        !   call mrsl(m1,prt(:),temp(:),rvls(:))
        !   do k = k2-1,m1
        !      vctr2(k) = c1
        !      vctr3(k) = c2
        !      vctr4(k) = c3
        !   enddo
        !   ki = m1 + 1
        !
        !   !	  if any ice phase microphysics are activated ....
        !   ! **(JP)** fatora expressoes logicas para fora dos lacos
        !   if (log3) then
        !   ! **(JP)** fim modificacao
        !      !  find level of -20 c.  assume ice saturation above this
        !      !   level.
	!
        !      do k = k1,m1
        !	  if (temp(k) .le. 253.16) then
        !	     ki = k
        !	     go to 10
        !	  endif
        !      enddo
        !      ki = m1 + 1
!10            continue
        !      call mrsi(m1-ki+1,prt(ki),temp(ki),rvls(ki))
!-srf 19/03/2005
!bug: Subscript out of range for array prt
!    subscript=0, lower bound=1, upper bound=132, dimension=1
! tempc < 253 sobre parte da Antartica.
!              ki = max(ki,k2)                                    !solucao 1
!              call mrsi(1,prt(ki-1),temp(ki-1),rvlsi)
        !       if(ki > 1) call mrsi(1,prt(ki-1),temp(ki-1),rvlsi) ! solucao 2
!srf
        !
        !
        !      do k = ki,m1
        !         vctr2(k) = ci1
        !         vctr3(k) = ci2
        !         vctr4(k) = ci3
        !      enddo
        !   endif
        !
        !   ! **(JP)** fatora expressoes logicas para fora dos lacos
        !   if (log4) then
        !   ! **(JP)** fim modificacao
        !      do k = k1,m1
        !         rc(k) = rcp(k,i,j)
        !      enddo
        !   else
        !      do k = k1,m1
        !         rc(k) = max(rv(k,i,j) / rvls(k) - .999,0.)
        !      enddo
        !   endif
        !
        !endif

        do k = k2,m1-1
           vctr1(k) = g / ((zt(k+1) - zt(k-1)) * rtgt(i,j))
        enddo

        ! **(JP)** fatora expressoes logicas para fora dos lacos
        !if (log2) then
        ! **(JP)** fim modificacao
        !   do k = k2,m1-1
        !      if (rc(k) .gt. 0.) then
        !         rvii = rvls(k-1)
        !         if (k .eq. ki) rvii = rvlsi
        !         en2(k,i,j) = vctr1(k) * (  &
        !              (1. + vctr2(k) * rvls(k) / temp(k))  &
        !              / (1. + vctr3(k) * rvls(k) / temp(k) ** 2)  &
        !              * ((vctr11(k+1) - vctr11(k-1)) / vctr11(k)  &
        !              + vctr4(k) / temp(k) * (rvls(k+1) - rvii))  &
        !              - (rtp(k+1,i,j) - rtp(k-1,i,j)))
        !      else
        !         en2(k,i,j) = vctr1(k)*((vctr12(k+1)-vctr12(k-1))  &
        !              / vctr12(k) - (vctr32(k+1) - vctr32(k-1)))
        !      endif
        ! 
        !   enddo
        !else
           do k = k2,m1-1
              en2(k,i,j) = vctr1(k) * ((vctr12(k+1)-vctr12(k-1))  &
                   / vctr12(k)                              )
        !          / vctr12(k) - (vctr32(k+1) - vctr32(k-1)))
           enddo
        !endif
        ! **(JP)** remove para fora do laco, permitindo vetorizacao
        !en2(k1,i,j) = en2(k2,i,j)
        !en2(nzp,i,j)=en2(nz,i,j)
        ! **(JP)** fim da modificacao

     enddo
  enddo

  ! **(JP)** removido de dentro do laco
  do j = ja,jz
     do i = ia,iz
        en2(lpw(i,j)-1,i,j)=en2(lpw(i,j),i,j)
     end do
  end do
  do j = ja,jz
     do i = ia,iz
        en2(nzp,i,j)=en2(nz,i,j)
     end do
  end do
  ! **(JP)** fim da modificacao

  return
end subroutine bruvais_OLAM

!     *****************************************************************


